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0.  Introductory  remarks 

Design  of  vibrating  structures  has  been  carried  out  on  trial 
and  error  basis.  No  comprehensive  theory  exists  at  the  present 
time,  although  some  progress  in  optimizing  the  design  of  a  single 
element  has  been  made.  P rage r  and  his  collaborators  worked  pri¬ 
marily  on  optimization  of  the  design  of  a  single  beam,  making  use 
of  the  Betti-Castigliano  deflection  formula  and  deriving  some 
theoretical  results  for  optimization  of  weight  for  a  given  de¬ 
flection  at  a  known  point  along  the  length  of  the  beam,  and  re¬ 
lated  problems  (see  [7]  or  [8]).  A  more  difficult  problem  of 
simultaneous  optimal  design  and  optimal  control  was  considered  by 
Komkov  and  Coleman  in  [4] ,  however  the  applications  of  their  work 
remain  limited,  and  no  general  optimization  theory  was  derived. 

A  different  static  optimization  was  pursued  by  Haug  [3], 

Haug  et  al  [2], .and  others,  who  attempted  to  derive  algorithms 
for  gradual  improvement  in  design  by  an  iterative  process.  So 
far  the  numerical  results  obtained  at  the  University  of  Iowa  are 
promising,  and  have  been  applied  to  some  problems  of  weapon  design. 

The  purpose  of  this  work  is  to  generalize  the  results  of  [3] 
via  sensitivity  analysis.  In  particular  to  derive  sensitivity 
(to  design  changes)  for  beams,  plates  and  general  structures, 
and  to  translate  these  theoretical  results  into  a  computational 
algorithm  using  some  form  of  the  steepest  descent  algorithm. 

A  particular  application  of  the  technique  derived  by  Komkov 
and  Coleman  ([43 )  is  given  to  the  problem  of  a  gun  tube  design. 
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1.  A  class  of  optimal  design  problems  described  by  a  distributed 
p ar ameter  system . 

k 

Consider  a  design  problem  postulated  on  a  subset  ft  £  R  ,  with 
a  local  coordinate  system  x  =  (x^ ,x2  *  • • • /X^) .  The  state  of  the 
system  is  determined  by  an  n-dimensional  vector  function  z  (x)  , 
which  is  an  element  of  a  Hilbert  space 

The  design  of  the  system  is  described  by  an  m-dimensional 
vector  function  u(x).  There  are  constraints  on  stress,  deflection, 
etc.,  and  physical  constrains  on  the  design  vector  u(x).  The 
set  of  m-dimensional  vector  functions  satisfying  the  constraints 
will  be  called  the  set  of  admissible  designs  and  will  be  denoted 
by  U  (x)  . 

The  state  vector  z  satisfies  a  system  of  differential  equations 

L  (u)  z  =  Q  (x , u)  ,  x  e  ft,  (1) 

where  the  forcing  function  Q  is  an  element  of  a  Hilbert  space  H^. 
Thus  L:  H^  -*  H^.  The  state  variable  z  also  satisfies  a  set  of 
boundary  conditions 

B  z  (x)  =  q(x),  x  e  3ft,  (2) 

Introduce  also  Hilbert  Spaces  H2,H4  of  functions  whose  domain 
is  3ft,  and  such  that  the  boundary  condition  (2)  is  given  as  a 
mapping  B:  H2  ■+  H^. 

While  the  same  symbol  z  is  used  in  relations  (1)  and  (2)  it 
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actually  depicts  different  classes  of  functions.  However,  z  re¬ 
garded  as  an  element  of  ©  H2  is  continuous  on  fl  y  3S2  and  no 
problems  will  arise  due  to  the  simplified  notation. 

The  system  (1)  and  (2)  will  be  referred  as  B.V.P.  (boundary- 
value  problem). 

Let  A  denote  the  operator 


H1  ®  H2  H3  ®  H4  • 


The  existence  of  a  unique  solution  of  B.V.P.  is  equivalent  to  the 
statement  that  A- ^  is  defined.  Well  posedness  of  the  problem,  in 
the  sense  of  Hadamard,  is  equivalent  to  the  statement  that  A  1  is 
defined  and  is  a  bounded  operator  from  into  .  This 

is  a  well  known  property  of  stable  structural  designs.  In  an 
important  class  of  structural  designs,  the  differential  operator 
L  also  turns  out  to  be  positive  definite  and  bounded  below 
(away  from  zero),  and  L-1  is  a  completely  continuous  operator. 

Also  of  interest  in  this  paper  are  designs  of  structural  mem¬ 
bers  subjected  to  buckling  and  to  natural  frequency  constraints. 

In  each  of  these  cases,  the  physical  problem  leads  to  an  eigen 
value  system 


L'  (u)  y  (x)  .=  ■  CM(u)  y(x)  ,  x  e  fl,  (3) 

BV  y  (x)  =0  ,  x  e  3 G,  (4) 

where  L'  and  M  are  symmetric,  positive,  bounded  below  differential 
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operators,  ?  is  a  real  eigenvalue,  and  y(x)  is  the  corresponding 
eigenfunction.  The  eigenfunctions  y(x)  are  normalized  by  the  con¬ 
dition 

<My,y>  =  HI  y HI  2  -  1, 

where  |||  •  1 1|  denotes  the  energy  norm  induced  by  the  inner  product 
JJX,Y]  =  <MX,Y>.  The  inner  product  <X,Y>  is  defined  in  the  usual 
way,  i.e.  <X,Y>  =  /^X(x)  Y(x)  dx. 

One  wishes  to  minimize,  or  approximately  minimize,  a  func¬ 
tional 

J  :  {H1  ©  H2  x  U}  -»■  R. 

Where  J  is  usually  either  an  L2  norm  of  Y,  an  energy  norm,  or 
a  possibly  non-linear,  but  positive  functional.  Analytically, 

J  is  of  the  form 

J(c;,  z(u),u)  =  hn(£)  +  /  gn  (z,x)dx  +  /  f  (z  (u)  ,u,x)  dx.  (5) 

The  design  problem  may  thus  be  reduced  to  selection  of  an  admissi¬ 
ble  design  u(x)  e  U  to  minimize  (or  approximately  minimize)  J, 
subject  to  the  constraints  (1)  -  (4)  and  performance  constraints 
of  the  form 

=0,  a=l,2,...r1,  (6a) 

<0,  a=r ■*■+!, .  . .  ,r  (6b) 
V 


f 

i|U  =  +  /  g„(2/X,u)dx  +  /f  (z,x,u)dx  ( 

■  a  8Q  a  ~  -  -  ~  12  a  ~  ~  ~  ~  ] 


cf>g(x,u)  ( 


A 


=0,  .3=1,2 , . . .  ?  , 


<0/  3=5+1,.../?. 


x  e  ft 


(7a) 

(7b) 


2.  Comments  based  on  computing  experience. 

The  pointwise  constraints  (7a)  and  (7b)  are  allowed  to  depend 
only  on  the  design  variable  u.  If  the  state  variable  z  were  pre¬ 
sent,  computational  difficulties  would  arise  later  in  the  develop¬ 
ment. 

Constraints  of  the  form  ri(x,u,z)  £  0,  x  e  ft,  arise  often, 
however,  and  must  be  treated.  It  is  readily  seen  that  the  point- 
wise  constraint  n  (x)  <_  0  is  equivalent  to  the  functional  constraint 


(h  (x)  +  |  n  (x)  I )  dx  =  0.  (8) 


Similarly,  <}>  (x)  <  an  Vx  eft  is  equivalent  to 

/Q{(4>(x)  -  aQ)  +  1 <t>  (x)  -  aQ  |  }  dx  =  0 . 

3.  Computation  of  Gateaux  derivatives  associated  with  an  optimal 
design  problem. 

The  standard  definition  (see  [6] )  of  a  Gateaux  derivative 
of  a  function  f:  B1  •+  B2  in  the  direction  of  a  vector  h  is  em¬ 
ployed.  Here  B^,  B2  are  Banach  spaces,  h  e  B^.  If  for  a  suffi¬ 
ciently  small  value  of  the  real  parameter  t 


f(xQ+th)  -  f(xQ) 


t  L  +  r  (xn,  th)  , 
x0 


(9) 
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where 

r(xn#th) 

lim  -  =0  ( 9a) 

t+0  t 

and  where  L  is  a  linear  operator  L  :  B..  ■+  B0 ,  then  f  is  Gateaux 

x0  X0  1  z 

differentiable  at  xn  (eB,)#  and  the  linear  operator  L  is  called 

u  j.  xQ 

the  Gateaux  derivative  of  f  in  the  direction  of  h.  If 
r (xn  # th) 

lim  - - - -  =  0 

t+0  t 

uniformly  with  respect  to  all  h  e  on  all  bounded  subsets  of  B2, 

and  if  L  is  continuous#  then  f  is  called  Frechet  differentiable 
x0 

at  Xq. 

Note:  the  left  hand  side  of  equality  (9)  will  be  called  the 
Gateaux  variation  of  f  in  the  direction  h,  computed  at  x^. 

In  this  paper  B^  and  B2  are  Hilbert  spaces  and#  in  the  develop¬ 
ment  that  follows,  B^  will  be  identified  with  the  Sobolev  space  Hk, 
while  B2  will#  in  most  cases,  be  the  real  line. 

The  minimization  problem  consists  of  finding  a  minimum  of  a 
functional  J  which  in  many  problems  of  structural  analysis  is  a 
quadratic  functional.  L2  norm  of  the  state  variables  an  energy 
norm,  or  possibly  other  norms,  subject  to  a  constraint  imposed 
on  the  state  variable,  with  which  must  satisfy  the  appropriate 
differential  equations  of  mathematical  physics.  That  is,  the 
problem  is  formulated  as  minimization  problem  for  the  functional 
J(u,z,C,x)  subject  to  differential  equations  Lx  =  £z.  For  example 


a2  d2 

Lz  =  ?z  may  be  of  the  form  — ylEI (x)  — y]  +  ;z  =  0.  Then  the 

dx^  dx~ 

minimization  problem  may  be  restated  in  an  enlarged  space  as  that 
of  finding  the  minimum 


.2  ,2 

A^  [<EI  (x)  ■■  2  /  — +  <£z,z>]  +  AqJ. 

dx  dx  u 


where  Aq,  A^  are  Lagrangian  multipliers. 

Hence  the  constraint  equation  is  regarded  as  one  of  the 
Euler-Lagrange  equations  which, subject  to  some  smoothness  assump¬ 
tions  on  the  state  variable  z, is  a  necessary  condition  for  the 
solution  of  the  minimum  problem  in  the  enlarged  space. 

Consider  the  functional  J(u,z ,x,£)  defined  by  (5)  and  its 
Gateaux  variation  6J,  where  6h  =  th  (for  sufficiently  small  (t) 
and  6h  =  (Sz ,  6u,  6  C)  .  In  a  purely  formal  fashion  one  may  write 


6  J  = 


8h“ 


3? 


"3z 


6z^> 


3ft 


3  f  3  f 

+  4^>  +<5V-'  6£>  +  r(S>- 

~  U  U  .  .  06 

(10) 


where  <  • , • and  <•  ,  • 
functions  whose  domains 


are  the  appropriate  inner  products  of 
are  3ft  and  ft,  respectively,  and 


l|r(h)|| 

lim  - - —  =  0. 

INI-o  l|h|| 


The  problem  arises  of  justifying  this  formal  approach  by  showing 
differentiability  of  the  operator  L'  with  respect  to  the  eigen¬ 
value  There  are  some  obvious  difficulties,  since  L'  is  an  un¬ 

bounded  operator.  Instead,  the  fact  that  A  ■*",  defined  by  la,  is 
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a  compact  operator  in  common  problems  of  structural  analysis,  can 
be  employed.  Hence,  the  eigenvalues  of  A  \  which  are  C  \  are 
a  discrete  subset  of  the  real  line,  i.e.  the  spectrum  of  A  ^  con¬ 
sists  only  of  isolated  eigenvalues.  One  may  now  involk  the 
theorem  on  bounded  invertibility  (see  Kato  [9]  pp.  196-197)  and 
a  well  known  result  on  perturbation  of  spectrum  of  A  namely 
that  the  spectrum  of  A-1  changes  continuously  with  A  \  such  small 
changes  being  interpreted  in  the  usual  operator  norm  topology ([§])• 
For  example,  consider  a  beam  equation  of  the  form 

^  /El  I  /IT  ^-3-  y  +  q  (x)y  =  f  (x)  +  Ly  (11) 

[dx2  dx2/J 

l2[0,1] 

subject  to  appropriate  boundary  conditions.  Let  qn  - — ->-  q, 
f  e  L2[0,1].  Then  R(?n)  -*■  R(£),  where  R  is  the  resolvent,  and 

lim  £  =  and  %  does  not  lie  in  the  convex  hull  of  6  implies 

n+°° 

that  the  solutions  of  (11)  converge  in  L2[0,1]  to  the  solution  of 
the  equation 


2  \ 

/I I  [/IT  2  y  +  q  (x)y  =  f  (x)  +  £y. 


The  key  in  this  discussion  is  the  avoidance  of  singular  pertur¬ 


bations. 


To  illustrate  this  remark  consider  the  following  two  examples 
The  fourth  order  beam  operator 


d2  a2 


EI(x)  , 

dx  i 
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with  boundary  conditions  of  either  free  support  at  both  ends,  or 

cante levered  at  one  end,  implies  continuous  changes  in  the  eigen- 
¥ 

values  with  a  continuous  change  in  EI(x).  The  operator 


y (0)  =  y'  (0)  =  0 

with  / 

I  y (1)  =  y'  (1)  =  0, 


(14) 


on  the  other  hand,  will  cause  singular  perturbations. 

(Note:  This  is  exactly  the  example  of  singular  perturbation 

given  in  [9]  page  436,  example  1.20) . 

That  is,  the  eigenvalues  of  (13)  are  stable  and  of  (14)  are  not. 
Note:  The  discussion  of  stability  of  eigenvalues  can  be  traced 
to  Lord  Rayleigh  (see  (10] ,  page  300) . 

It  is  assumed  in  the  use  of  formula  (10)  that  A-1  is  compact 
and  that  perturbations  of  eigenvalues  of  A-1  (and  of  A)  are  stable. 
That  is,  one  is  dealing  with  a  common  case  of  stable  structural 
design. 


Manipulation  of  formula  (10) . 


The  meaning  of  <$£  and  its  relationship  to  <$u  needs  to  be 
determined  next  in  the  interpretation  of  formula  (10) .  Regarding 
SC  as  the  Gateaux  variation  of  the  first  eigenvalue  of  the  system 
L'y  =  CMy ,  and  assuming  that  L'  and  M  depend  holomorphically  on  u, 
it  is  possible  to  apply  Arons zajn's  theory  of  general  quadratic 
forms  and  justify  formal  differentiation  of  the  associated  Rayleigh 


quotient.  (see  Kato  [9],  VII  -  83,  pp.  421-422). 

The  first  eigenvalue  of  the  system  is  equal  to  the  value  of 
the  Rayleigh  quotient 


C-,  = 


<L’Yi-yi>a 
1  <%'  Xi>  ’ 

where  (x)  is  the  corresponding  eigenfunction.  Assuming  symmetry 
of  L*  and  M  and  positive  definiteness  of  M,  and  using  formal  rules 
of  Gateaux  or  Frechet  differentiation  (see  for  example  the  intro¬ 
ductory  notes  in  [11])  for  the  Rayleigh  functional,  one  obtains 

=  (MYl'Yl)‘a  {4L'Yl>h'Yl)  -  ^(“Yl'h'Yl)  +  (<L'Yl  ' 


cMyx) ,y1  \  >  / 


(15) 


where  subscript  h  denotes  Frechet  differentiation  with  respect 
to  a  typical  component  h  of  the  vector  u. 

Expressions  of  the  type  (L'y)h  have  to  be  carefully  inter¬ 
preted,  since  y^  and  the  map  L1  both  depend  on  h.  If  both  are 
Frechet  differentiable,  the  following  rule  is  easily  justified 
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or  equivalently 


=  LhY  +  Ly 


where  L^y  has  to  be  interpreted  as  the  Gateaux  variation  of  the 
operator  L  with  y  regarded  as  fixed,  while  y^  is  the  corresponding 
variation  of  y . 

Observe:  L (hQ+th) y (hQ+th) -L (hQ)y (hQ)  =  [L (hQ+th) y (hQ+th) -L (hQ ) y (hQ+th) ] 

+  [L (hQ)y (hQ+th) -L (hQ)y (hQ) ] . 

The  Gateaux  variation  of  the  first  eigenvalue  6  is  given  by 


Since  y^  is  the  eigenfunction  of  (3)  the  last  term  can  be  omitted, 
and 


Note 


(17) 

/9L'yl  \ 

5  \  1u" — '  Yl/  is  not  a  scalar,  it  is  a  vector  in  the 


space  dual  to  u.  The  same  remark  applies  to  the  vector 
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5.  An  adjoint  problem. 

Since  the  generalized  displacement  (or  state)  vector  z  is  an 
element  of  the  Hilbert  space  ©  H2,  its  dual  (the  generalized 
force)  is  an  element  of  the  same  Hilbert  space.  The  duality  is 
maintained  for  sake  of  physical  interpretation. 

Let  XJ  be  an  element  of  the  dual  space,  satisfying  the  differ 
ential  equation 


J  9f0 

L*X  =  —  in  > 


where  L*  is  a  formal  adjoint  of  L,  satisfying 
<W,Lz\  -  <L*W,z\  +  /cw,z  > 

V  n  x  ~  a  x  ~  an 

where  C  is  a  linear  operator  C:  H2  -*\  H2. 

\J  satisfies  boundary  condition 


for  every  admissible  6z,  which  is  certainly  satisfied  if 


CXJ  =  on  38  . 


/l«?,  xj\  .  ^y  +  <cxj. y 

\  '  o  "  3 


Hence 
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3ho  /  3gn  v  /  9fo  \  /  9fn  \ 

5J  -  rr  65  +  <  v*’’  6~/  +  \  ^  6?  /  +  <  5/'  *."/ 

an  77,;  ~7  a :  a 

=  TT  6C  +<vS\T'6~^  +^J(“5if“  "  ^  +  ^  Fz^'6?) 


-  /ex^zV  -  ■:,*.;:•■■ 

;  N  ~  an  ... 

._  3h0  /.  .  /  /3f0  ,  \J  / 9  (Lz)  3Qv\  x.;\ 
=  TT  u  +  \  (3u“  +  ~  (Tu“  ■“■  tt]) '  62/ 

..•■■/■■ '  ~  .  .  ;  r 

3h  .  ....  '  .  ,2" -■ 

“  jt-  6?  +  k  -•  &~y  • 


7  _  9f0  .  ,J  73  (Lz)  30s 

i  -  +  2  (-^ —  <Tu 


is  the  sensitivity  vector  if  ^  6C  2  0 ,  i.e.  the  first  natural 
frequency  does  not  change  significantly  with  the  admissible  changes 
in  the  design. 

3hn 

If  the  term  t—-  6?  can  not  be  ignored,  a  substitution  of  (17) 
into  (23)  yields 


an.  >  v  V  /  a  (L'y. ) 

4 J  -  TT  <  MXl  -Xi>  {\  1  5u  Xl  -  l 

*  { 5s)„  -  7  7  • 


3 (My. ) 


r,l»  5u)  "'}■ 

1  ~  V  Q  . 


A  is  the  sensitivity  vector  for  optimization  of  J. 

V  It  is  Significant  that  no  boundary  terms  (i.e.  products  in 
H2)  occur  in  the  formula  (25) . 
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6.  An  example.  Consider  the  dynamic  problem  of  thin  plate  theory 
represented  by  the  equation 

2 

V2  (D  (x,y)  V2W)  +  =  Q(x,y)e:|a)t,  x  e  O',  t  >  0  . 

dt 

After  separation  of  variables,  this  is  reduced  to  an  eigenvalue 
problem 


x 

(L*  +  CD  W(x)  =  Q (x)  x  =  •  (26) 

y 

U  J 

..  ,  - 

where  L'  =  V2(D(x,y)V2)  and  M  si  is  the  identity  operator. 

Design  a  plate  of  a  given  weight  (say  unity)  such  that  the 
average  kinetic  energy  takes  an  approximately  extremal  Value.  The 
dep)g)gn  parameter  u(x)  is  the  thickness  of  the  plate.  The  plate 
is  clamped  on  the  boundary.  In  fact  this  problem  is  equivalent 
to  looking  for  a  low  approximate  value  of  the  sensitivity  vector 
A. 

The  cost  functional  J  is  identified  as 


j  =  y  fn  u(x)  dx  +  kC 
6  J  M  ~ 


(27) 


where  k  is  a  known  constant.  (Recall  that  £  is  proportional  to 
w2) .  Identifying  h^(S)  as  k£,  gg  =  0,  fg  =  yQu,  the  sensitivity 
vector  is  computed  as 


A1  -  y. 


♦  k  <y,,Yl>  -1  {<V2  (2^-  V2yi  +  -ilp 

1  x/  N  3  (1-v2)  1  12 (1-v  ) 


9  (V  Yl) 


'yl) 


'  SYi 

Supposing  that  5^  and  y1(x)  are  known  and  is  computed,  the  sen 
sitivity  vector  (28)  can  be  computed  in  terms  of  known  physical 
quantities.  An  optimum  design  requires  near  zero  sensitivity. 

In  order  that  Uq  is  optimal  it  is  necessary  for  AJ  to  vanish. 
However  equation  (28)  is  difficult  to  solve  for  Uq  and  gradual 
improvement  approach  is  generally  used  for  suboptimal  design. 

If  AJ  is  not  zero  then  following  a  steepest  descent  appraoch  one 

'•  ■  \  j  -  .  ■ 

could  choose  6u  in  the  direction  of  -A  .  This  is  locally  optimal 
direction.  In  general  the  constraints  improved  by  physical 
consideration  may  not  be  satisfied  by  altered  design  vector  yQ  + 
<$u.  Various  versions  of  the  gradient-projection  technique  have 
been  used  to  correct  this  defect. 

For  an  arbitrarily  chosen  design,  the  decrease  in  the  mag¬ 
nitude  of  the  sensitivity  vector  indicates  the  degree  of  improve¬ 
ment  in  the  corresponding  step  of  the  steepest  descent  algorithm. 

7.  A  theoretical  foundation  of  the  gradient  projection  algorithm. 

Suppose  that  one  wishes  to  minimize  the  cost  functional 
J(z(u),u)  subject  to  constraints  of  the  form 

A  z (u)  =  b, 

where  z  belongs  to  a  Hilbert  space  H^  b  to  H2,  and  A:  H1  ^  H2 
is  a  bounded  operator  with  closed  range,  whose  domain  is  dense 
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in  H^.  For  subsequent  analysis,  N (A)  denotes  the  null  space  of 
A  and  u  denotes  an  m- dimensional  vector,  called  the  design. 

It  is  assumed  that  J(z(u),u)  =  J(u)  is  a  Frechet  differentiable 

A 

functional.  J* (u)  will  denote  the  gradient  of  J.  In  the  pre¬ 
vious  discussion  of  sensitivity  (parts  1  -  5)  of  this  paper  J' (u) 
was  identified  with  the  vector  AJ. 

Suppose  an  initial  choice  of  design  u  is  made,  uQ,  with  the 

a  ; 

corresponding  value  J(Ug)  computed.  Az(Ug)  =  b(Ug)  is  (approxi¬ 
mately)  satisfied.  Choose  <Su,  so  that  |[|6u ||p  £  e  (in  fact  the 
choice  will  be  ||(5u. |||  =  e,  unless  a  change  or  a  refinement  of  step 
size  is  required).  If  no  constraints  were  present,  the  steepest 

A 

descent  condition  would  be  (for  some  constant  k  >  0)  6u  =  -kJ' (un) . 

A 

Note  that  J' (u)  is  an  m-dimensional  vector  in  the  space  dual  to  U. 
Mathematically  these  spaces  are  indistinguishable,  i.e.  isometrically 
isomorphic,  even  though  the  physical  dimensions  are  different. 

For  this  reason  the  symbolism  6u  =  -k  J' (u)  makes  sense.  Of 
course  in  the  general  case  such  choice  of  change  in  the  design 
u  will  violate  the  constraints. 

For  the  sake  of  convenience  the  constraint  can  be  put  in  the 
form  of  a  functional  equality 

tp  (u)  =  0,  ip:  U  -*■  R. 

For  exanple,  A  z  (u)  -  b(u)  =  0  can  be  replaced  by  *Mu)  =  ||A  z(u)  - 
b  (u)  ||  =  0. 

In  discussion  of  this  technique  the  following  problem  will 


17 

be  considered  first:  The  functional  J(u)  is  to  be  minimized, 
subject  to  a  constraint  Au  =  0,  where  A  is  linear  operator, 

A:  Hx  -»  H2. 

In  this  case  the  algorithm  describes  a  consecutive  steepest 
descent  and  projection  on  a  constraint  surface  in  a  single  step. 

A 

J' (u)  can  be  uniquely  represented  as  a  sum  of  a  vector  in  N(A) 
and  a  vector  in  the  range  of  A*;  (Recall  that  R(A*)  =  N  (A) 1. ) 

-J' (u0)  =  Y0  +  A*no  (7.1) 

where 

Ay0  =  0.  (7.2) 

Hence , 

0  =  -AJ' (un)  +  AA*n  (7.3) 

Since  AA*  is  invertible  on  (N  (A) ) 1  (and  Rq  e  (N  (A) ) 1) ,  rq 
can  be  computed  as 

ri0  =  -  (AA*)  "*1  A(J'  (uQ)  )  .  (7.4) 

Hq  is  the  orthogonal  projection  of  -J' (u)  on  the  constraint  hyper- 
surface  Au  =  0. 

Setting  u^  -  Uq  +  ay^  completes  the  first  iterative  step, 
where  a  +  ay^  is  a  constant  chosen  to  satisfy  constraint  on  step 
size,  0  <  a  £  1.  The  iteration  proceeds  according  to  the  formulas 

Y  =  -[I  -  A*(AA*)“1A)  J'(u  )  (7.5a) 

n 
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u  J_,  =  u 
~n+l  ~n 


(7.5b) 


If  the  problem  has  a  unique  solution,  this  iterative  technique 
converges  to  that  solution  (See  [12]). 

Critique.  In  most  structural  design  problems,  operator  A  is  not 
linear.  A*  is  generally  undefined,  and  the  above  technique  appears 
to  be  difficult  to  apply.  However  if  the  constraint  functional  is 
Frechet  differentiable,  then  locally  (for  small  variations  of  the 
parameter  vector  u)  the  constraint  can  be  approximated  by  a  linear 
one.  That  is  the  hypersurface  \ p  =  0  can  be  locally  approximated 
by  a  supporting  hyperplane.  However  the  steepest  descent  step 
is  projected  on  such  a  hyperplane  and  not  on  the  constraint  sur¬ 
face,  which  means  that  the  constraint  is  not  satisfied,  but  is 
only  approximately  satisfied. 

An  additional  correction  could  be  applied  at  each  step  (or 
only  at  certain  steps)  of  this  iteration  to  insure  exact  complience 
with  the  constraints.  For  example,  a  step  6u  =  -8^' ,  0  <  | 3 | .<  1 

A 

can  be  put  in,  with  J  3 1  estimated  to  be  ||i|>||  •  ||t£’  \\  ,  to  reduce 

close  to  zero  the  value  of  ||^|j.  8  is  positive  if  ip  is  positive, 

eind  negative  if  ip  is  negative. 

Remark :  In  most  engineering  problems  exact  complience  with  con¬ 

straints  is  not  very  important,  since  the  constraint  conditions 
are  inaccurate  and  were  derived  following  some  Simplifying  assump¬ 


tions 
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7—.  Haug's  variant  of  the  steepest  descent  algorithm  for  structural 
design  optimization 

In  [3]  Haug  has  employed  a  technique  which  can  be  regarded  as 
a  variant  of  the  gradient  projection  method.  Further  refinements 
and  some  changes  in  this  approach  were  given  by  Haug,  Arora  and 
Matsui  in  [14] .  The  problem  involves  minimization  of  a  functional 
J  =  ^0(u'z'£)'  subject  to  inequality  constraints  <£  (u)  <  0,  i|k  (u)  <_  0 

j  =  1/2,.'. .  ,k;  as  in  (8)  pointwise  constraints  being  replaced  by 
functional  constraints.  As  usual 

5  4>  =  ,  Sip  m  6u>  +  (ipz,  Sz)  +  (7.6) 

Explicit  dependence  of  on  6z  and  6c  is  eliminated  by  perturbing 
the  state  equations  (1)  and  the  corresponding  Hamiltonian  functional 
FL (a,u, X)  =  <L(u)z,  X>  =  <Q  (x,u)  ,  X>  =  j|)(u,z)  ,  (7.7) 

where  X  is  physically  a  dual  element  of  Q (x ,u) .  Since  the  dual  of 
H3  is  H3,  X  is  an  element  of  the  Hilbert  space  H3.  Taking  a  Fr^chet 

derivative  of  FL  with  respect  to  X  one  recovers  the  equation  (1). 
Differentiating  with  respect  to  z,  one  arrives  at  the  adjoint  system 


L*  (u)  X  =  -ipz  (7.8) 

This  leads  to  identities 

A  ,  L6  z  )  +  (Lz )  6u  )  =  0  (7.9) 

~  ~  ~  U  ~ 

( L*X ,  6z  )  +  CX,  (Lz)  <5u>  =  0  (7.10) 

.  ~(4>z,  <5z  >  +  <X,  (Lz)u  <5u)  =  0  (7.11) 


Since  (Lz)  ;6u  =  -L  6z,  dependence  of  Sip  on  Sz  is  eliminated.  If 
y  is  an  eigenfunction  of  the  system 
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L(u)y  -  ?  M(u)y  ,  (7.12) 

where  as  before  M  is  positive,  definite,  and  symetric,  one  derives 
an  identity 

(Y,  ( (Ly)u  -  C(My)u)  6u>  =  Since  <y,M(u)y  >  =  1 

6\p  =  <Su>  +  <  X,  (Lz )  u  6u>  +  \p^  <  y,  ( (Ly)u  -  C  (My  )u  6u> 

=  (  X  ,  <$u )  £  Aty  .  (7.13) 

Here  one  is  approximating  the  non-linear  terms  by  first  order  lin¬ 
ear  approximations. 

A  constraint  on  the  magnitude  of  6u  is  imposed  by  requiring 
2 

<W6u,  6u)£  §  where  the  operator  W  (the  weighing  operator)  is  posi 

2 

tive  definite,  and  symetric.  §  is  chosen  sufficiently  small  for 
the  particular  problem. 

Application  of  the  multiplier  rule  leads  to  the  set  of  equations 
X°  +  y  •  <J>  +  y*X  +  2v*W  <5u  =0 

d)  •  6u  -  A  4>  —  0  (7.14) 

/  •  {< A ,  6u>  - Aty}  =  0  (7.15) 

v  {<  Wfiu,  <5u>  -  $2}  =  0  (7.16) 

n 

Here  (  .  )  means  the  product  A  •  B  =  £  A.  B.  as  used  in 

..  i=l  1  . 

the  system  of  equations  (7.14)  -  (7.16),  The  equations  (7.7)  - 
(7.11k  can  be  formally  solved  for  6u,  while  (7.14)  -  (7.16)  can 
be  used  to  eliminate  the  multipliers. 

The  following  formulas  are  derived. 


6u2  -  W_1  LI  -  1  S'.w'1] 

,x[*  “*}  {**  •  Mk)J  +  W_1*u  (  *uW_1*u)  A*' 

where 

M»j  =C'.v‘lf  I  -♦a(*u»'1*u)  1  «U«_1J  A°), 

M«  -(i''w'1*u(*uw‘1+u)  1  *♦>  ' 

and 

M«  =  (*'  W_1  t1  -  *u(*uW"^u)  1  *uw_1]  *)■• 

“jj  -  <A°,  Vf1  [  I  -  ^  f  +UW-1^  )  1  »0W_1J  A0)  } 

while  fill  is  given  by 

6u  =  6u^  +  6u^, 

2  i 

6u  is  basically  the  constraint  correction,  6 u  gives  a  step  in 

the  direction  of  negative  gradient  projection.  Simultaneous  ap- 
1  2 

plication  of  6u  ,  <5u  constitutes  Haug's  version  of  the  modified 
gradient  projection  method. 

Applications  to  truss,  beam,  and  plate  design  can  be  found 
in  [14] . 


Optimal  design  of  statically  determinate  structures. 


In  the  statically  determinate  case,  the  loads  and  moments 
transmitted  by  the  members  are  independent  of  the  geometric 
design  of  each  member.  Hence  if  the  design  vector  influences 
only  the  cross-section  property  of  each  member,  the  sensitivity 
vector  is  relatively  easy  to  compute. 
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As  an  example,  consider  a  statically  determinate  beam  in  pure 
bending.  The  deflection  is  given  by  the  Betti  formula 

»<*>  -  /  wrfrr-  a«  =  *  <“*”»  •  (8-1) 

o 

where  M  (the  moment  applied)  and  m  (the  moment  at  x  due  to  a  unit 

load  at  £)  are  independent  of  the  design  u(£).  Using  rules  of 

3W 

Frechet  differentiation  the  sensitivity  operator  is  computed  as 

If  -s  ((MI-2)  *  m>  ||  (8.2) 

In  the  statically  determinate  case  it  is  therefore  unnecessary  to 
eliminate  since  all  terms  in  the  formula  (8.2)  above  are 

known. 

The  sensitivity  vector  AJ  can  be  identified  with  the  Freshet 
derivative 

i  (MI-2  (u)  *  m)  2|M,  , 

'  i 

if  no  constraints  are  imposed.  A  trivial  result  is  obtained  by 
setting  AJ  =  0. 

If  a  weight  constraint,  and  or  minimum  cross-section,  or 
maximum  stress  constraint  is  imposed,  the  optimum  design  solution 
becomes  non-trivial. 
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In  case  where  a  functional  relation  ij>(u)  can  be  derived  for 
the  constraint  i p  <_  c,  the  necessary  condition  for  optimization  of 
the  design  u  is  of  the  form 

sr  (MI'2  <“>  *>>  +ii5i  =  9- 

This  remark  is  illustrated  by  the  following  design  problem. 

9.  A  design  problem  for  an  optimal  shape  of  a  gun  barrel. 

The  gun  tube  design  can  be  regarded  as  a  one  parameter 
family  of  maps  h(x)  e  '■+  max  1 0  (x)  j  e  R.  is  a  class  of 

functions  obeying  certain  constraints,  and  it  will  be  called 
the  Class  of  admissible  thickness  design  functions.  Physically 
h(x),  x  e  [0, A],  represents  the  thickness  of  the  gun  tube.  Rq 
will  denote  the  caliber.  . 
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We  formulate  the  optimization  problem  as  follows; 

Find  an  function  ft(x),  x  e  10, &.]  such  that  min  (max 

|y«  $(X)  ,x)  | )  £  min  (max  |  y '  (h(x)  ,x)  | )  for  any  other  choice  of 

h (x)  e  L9 [0, &] ,  with  the  constraint  h(x)  £  hQ  tfx  e  10, i) ,  the  max- 

M 

imum  stress  constraint  amax  £  a0/  an<^  constraint  W  £  WQ. 

The  gun  tube  is  assumed  to  be  a  cantilevered  beam  subjected  to 
the  a  priori  given  accelerations  in  the  direction  of  its  axis 
and  in  the  plane  perpendicular  to  its  axis.  If  more  than  one 
optimal  design  exists,  then  the  design  which  is  optimal  according 
to  our  criterion  (i.e.  it  minimizes  max  dy  *  and  which 

also  results  in  smallest  total  weight  is  deemed  the  optimal  design. 


Mathematical  formulation  of  the  problem 

Since  the  beam  is  simply  cantilevered,  and  we  have  a  stat¬ 
ically  determinate  case,  we  can  use  Prager ' s  approach.  The 
bending  moment  M(x)  and  the  moment  m(x-xQ)  (at  x  due  to  a  unit 
load  positioned  at  xQ)  are  independent  of  the  design.  Castigliano' 

'■  .  '  jj  • 

theorem  can  be  rewritten  in  the  form: 


y(xft)  =  E”1  /  (I(x))"xM(x)m(x-x0)dx 


-1. 


(9.1) 


We  seek  to  minimize  max|h' (xQ) |  -  max  {  /(I(x)  M(x)  m' (x-xQ)dx}, 


xQ  e  [0 


,i],  subject  to  constraints  /  p  dx  £ 


W, 


(9.2) 


h(x)  >  hn. 


p  amax 

JL 

M(X)  •  C  (X) 

A  (XT- 

T 

I(x) 

<  a  , 
—  o  ' 


(9.3) 


2  2 

A  (x)  =  7r(R0  +  h(x) )  -  ttRq 


(9 .  4a) 


where 
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4  4 

7Tp(Rn  +  h(x) )  pTtRn  h 

I  (x)  =  - - - -  -  - (9.4D) 

4  4 

p  is  the  material  density, ■ ama  :  is  maximum  acceleration  exper¬ 
ienced  by  the  gun  mount. 

The  maximum  force  (per  unit  length)  due  to  inertial  load  is 
q(x) ...  ■  -amax1'p(h  +  R0)2  -  R0)'  (9,5) 

Note: 

In  this  computation  we  shall  ignore  the  possibility  of 

Tp  Mp 

torsional  vibration  and  stresses  and  use  the  formula  a  =  +  tt— . 

It  is  done  for  reasons  of  simplicity.  However ,  incorporating 
the  effects  of  torsion  stress  into  this  computation  is  quite 
straightforward. 

The  static  deflection  is  given  by  the  formula: 

a 

y  (Xq)  =  E  1  /  (lR0+h(x)  ]4-Rq)  }-1*M(x)  *m(x-x0)dx.  (9.1a) 

It  can  be  shown  that : 

clx  x  _Xq  “  E_1  / '  ( [RQ+h]  4-Rq}-1M(x)  (m(x-xQ))  dx  (9.6) 

Let  us  denote  the  quantity  (9.6) .by  9 (Xq) .  The  problem  is  now 

reduced  to  the  finding  of  necessary  conditions  for  minimization 

of  max  |9(Xq)|  subject  to  constraints  (9.3,). 
x0 

The  formulas  for  M(x)  and  m(£) (£=x-xn)  are  given  respectively 

by:  M (x)  =  /  q(?) •  (£-x)d5  (9.5a) 

x 
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where  q(x)  is  given  by  equation  (9.5) 

m(5)  =  u(-5>#  where  U(5)  is  the  unit  step  function 

U(5)  =  0  if  K  <  0 

=  1  if  K  >  0. 

U«  (5)  =  6(5)  (6  will  denote  the  Dirac  delta  function) . 

Hence  we  have 

0  (x)  -  E_i  (M(x) KRo+h)  4-^]  )  "1*  <$(*) 

=  (I^e.)"1  M(x)  •  ([RQ+h(x)]4  -  Rq)  •  (9.6a) 

(*  denotes  the  operation  of  convolution) . 

We  consider  two  possibilities: 

(a)  max  1 0  (x)  |  occurs  at  an  interior  point  of  I0,£].' 

x  e  [0,&] 

„(b)  |  0  (fc)  |  =  max  |  0  (x)  |  . 

xe I0,£] 

The  case  (b)  will  be  ignored  for  the  time  being,  since  6(1)  is 
easily  computed  (for  a  given  design)  and  Prager's  arguments  ([3]) 
can  be  duplicated  with  only  minor  modifications.  The  case  (a) 
is  more  complex. 

Our  discussion  will  concentrate  on  minimization  criteria 

for  the  case  (a),  and  only  a  final  comparison  shall  be  made  to 

check  if  I  0  U)  I  exceeds  max  |  0  (x)  |  ,  b  <  &. 

xe I 0 ,b ] 
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A  necessary  criterion  for  optimization  of  an  interior  iniaximum 
(case  (b) )  . 

Following  the  Komkov- Coleman  .technique  outlined  in  14]  we 
equate  to  zero  the  Frechet  derivative  of  the  functional: 


<Mh)  =  1.0(h)  +  X-^h-hg)  +  X  2  {--vrcr-  + 


I M  (h)  I  c  (h) 


+  X  3  3-  (0(h))  +  X.(/  A  (h  (x)  )  dx  -  W  ). 


(9.7) 


$h  =  X0(^|2-)"1{(  (x*(-27tph)  )  [  (R0+h)  4-R^]-(7rpx*  ((h+RQ)2-R2) ) 


4  (R+h)  3  +  X  i  +  X4Ah  +  X2  (sign  (M)  •  (f%h)  +  X3  ^  (@h)  =  0  (9.8) 


The  Frechet  derivative  of  0(h(x),x)  denoted  by  0^  is  exactly 
the  first  term  of  the  above  expression,  i.e. 


6h(x)  =  (£|P)-1{((x*(-2Trph))- [(RQ+h)4  -  R4]  +  (-7TPx* 


((h  +  RQ)2  -  Rq) * 4 (RQ+h) 3) . 


(9.8?) 


(as  before  *  denotes  the  operation  of  convolution) .  X4,X1,X2 

obey  the  inequalities : 


X  j  £  0,  X2  <  0,  X„  4  0 « 


(9.5) 


XQ,  X^,  Xji  X4,  are  constants,  but  X^  must  be  regarded  as  an 
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unknown  function.  The  relation  (9.8)  and  inequalities  (9.9)  deter¬ 
mine  a  necessary  condition  for  a  maximum  of  0  (x)  subject  to  the 
constraints  (9.3).  ’’ 

The  necessary  condition  for  optimization  of  max|9(x) |  (x-an 
interior  point  of  [Q,Jt))  given  by  the  relation  (9.8)  states  that 
(eh)  lies  in  the  forward  cone  of  the  three-dimensional  sub¬ 
space  of  I<2  [0/£]  spanned  by  @h,  Ah,  and  by 

{<si*nM<h)Wf£>h  -  V'  (?-10) 

1  h  A2  (h)  h 

where  subscripts  denote  Frechet  differentiation.  ^^(0^)  "*‘S 
computed  directly . 

^_9h  =  (I^E)"1  {  [-2irph(R0+h)  4-Rq]  +  [  (x*  (-2?rph) )  •  ( (RQ+h)  4-R4] 

-(47rp*((h  +  Rq)2-R2))  *  (Rq  +  h)3]  -  3  (TTPx*(h+R0)2-R2)  ? 

(RQ  +  h)2  h’ (x) }  .  (9.11) 

We  observe  that  ©h(x)  satisfies  a  first  order  differential  equation 
(8) ,  which  is; 

Vh  +  x3  a;  <V  =  "xi^2  (sign  M(h)  ■(^h)"x4V: '  <9-12> 

X3  =  0  corresponds  to  the  case  when  maximum  of  0 (x)  occurs  at 
the  boundary  point  x=£. 
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Moreover,  A 3(x)  =  0  for  some  x  e  |Q,A)  implies  that  x  is  a 
singular  point  of  the  equation  (8),  and; 


-A 


eh(x)  = 


X,  M(h)  (Rn+h) 

t—  (sign  M(h(x) )  •  (- - -), 

*•  Kh)  1 


-  T--A, 
x— x  *0  h 


x=x 
(9.13) 


However,  6^  is  a  known  function  of  h(x)  (see  equation  8a)  and 
the  numbers  obtained  by  computation  (9.13)  and  (9.8a)  after  substi¬ 
tution  of  some  design  function  h(x)  will  generally  fail  to  be 
equal.  The  possibility  that  they  could  be  equal  will  however 
remain,  and  an  additional  check  will  be  added  to  the  computational 
algorithm  to  take  care  of  that  possibility. 

The  preceeding  example  could  be  rewritten  to  incorporate  the 
results  of-  sections  4  and  5.  Assuming  that  the  gun  barrell 
vibrates  with  a  natural  frequency  the  discussion  above  can 

be  modified  by  incorporating  the  additional  term  similar  to  the 

2  '  '  d2 

second  term  in  equation  (28)  ,  where  is  replaced  by  — *•,  and 

dx*2 

Y1  by  (x)  which  is  corresponding  eigenshape. 

10.  An  example  of  computation 

Using  the  theoretical  approach  of  sections  4  and  5.  E.T.  Haug 
S.S.  Arora  and  K.  Matsui  have  programmed  an  optimization  algorithm 
for  a  simply  supported  beam  of  rectangular  cross-section  with 
width/height  ratio  of  1:2,  and  I  (x)  =  ^-u  (x)  where  u(x)  is  the 
design  parameter. 

The  loading  is  shown  below 
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1100  in  lb 


Figure  2 : 


The  lower  bound  on  cross-section  was  selected  to  be  .3645  in  2. 

The  displacement  constraint  is  A,  the  first  natural  frequency  is 
denoted  by  Cq.  The  cost  functional  is  the  total  volume  of  ma- 
terial.  (Here  E  =  10  psi,  p  =  .00208  slugs/in  .  The  computation 
was  performed  on  IBM  360  at  the  University  of  Iowa.  Convergence 
was  rapid-no  more  than  35  iterations  were  required  in  any  com¬ 
putation.  The  optimum  design  is  illustrated  below  for  various 
values  of  A, 

It  is  interesting  to  observe  how  complexity  of  design  is 
influenced  by  the  choice  of  AQ  in  the  constraint  A  £  AQ. 

A  similar  phenomenon  is  observed  if  constraint  on  the  value 
of  maximum  stress  is  imposed,  instead  of  the  maximal  displacement. 
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11.  Conclusions . 

The  technique  derived  in  sections  4  -  6,  and  the  more  spec¬ 
ialized  results  of  sections  8  and  9  have  been  shown  to  be  adaptable 
to  numerical  optimization  schemes  which  compare  favorably  with 
other  methods,  such  as  adaptation  of  Pontryagin's  maximality 
principle  proposed  in  13 J.  Direct  solutions  of  the  design  problem 
with  complex,  non-linear  constraints  represents  a  new  result  in 
this  area  of  engineering  analysis.  Even  the  optimization  problem 
discussed  in  section  8,  and  illustrated  in  example  9  can  not  be 
derived  by  a  direct  application  of  Prager's  result.  The  problem 
discussed  there  is  not  the  dual  of  Prager's  problem  17],  nor  can 
it  be  deduced  as  a  consequence  of  the  works  of  Prager's  collaborators 
(which  have  not  been  quoted  in  the  ref erences) to  best  of  our 
knowledge . 

Computing  experience  with  the  approach  suggested  in  this 
paper  have  been  satisfactory.  In  fact  satisfactory  numerical 
results  have  been  obtained  in  cases  where  theoretically  the 
technique  derived  here  is  not  justifiable. 

Additional  theoretical  work  is  clearly  needed.  Perhaps,  a 
more  sophisticated  approach  based  on  abstract  variational  prin¬ 
ciples  (see  [13])  may  by-pass  some  of  the  difficulties  encountered 
in  justifying  purturbation  type  arguments,  or  arguments  based  on 
vanishing  of  Frechet  derivative  and  the  use  of  Lagrangian  multiplier 
rule  in  problems  of  structural  mechanics. 
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